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ABSTRACT 

This  thesis  describes  a  Matlab  computer  program  that  implements  Amott’s  fonnulation  of 
thennoacoustics  (W.P.  Amott,  et.  al.,  "General  formulation  of  thermoacoustics  for  stacks 
having  arbitrarily  shaped  pore  cross  sections,”  J.  Acoustic.  Soc.  Am.,  90(6),  3228-3237 
(1991)].  The  program  calculates  the  resonance  frequency  and  the  quality  factor  of  a 
thermoacoustic  prime  mover  below  onset  of  self-oscilation.  The  results  of  this  analysis  are 
compared  to  measured  values  for  both  a  closed  end  and  an  open  end  prime  mover  and  to 
predictions  of  a  standing  wave  aiudysis  of  prime  movers  (A.A.  Atchley,  "Standing  wave 
analysis  of  a  thermoacoustic  prime  mover  below  onset  of  self-oscillation,"  J.  Acoustic.  Soc. 
Am.  ,  22(5).  2907-2914  (1992)]. 
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I.  XMTRODOCTIOH 


Thermoacoustic  heat  transport  is  a  process  through  which 
an  acoustic  field  generates,  or  inversely  is  generated  from, 
a  flow  of  heat.  Like  every  engine,  a  thermoacoustic  engine  can 
be  configured  as  either  type  of  classical  heat  engine,  a  heat 
pump (or  refrigerator)  or  a  prime  mover.  The  purpose  of  this 
thesis  is  the  computer  implementation  of  a  general  formulation 
of  thermoacoustics  developed  by  Arnott  et  al  (Ref.  1]  using 
Matlab.  In  particular,  we  find  the  theoretical  values  of  the 
quality  factor  and  the  resonance  frequency  of  a  prime  mover  as 
a  function  of  the  temperature  difference  applied  across  the 
prime  mover  stack.  The  program  determines  the  impedance  and 
normalized  pressure  distribution  along  the  prime  mover 
beginning  from  a  known  value  of  specific  acoustic  impedance 
and  normalized  pressure  amplitude  at  one  end.  It  finds  the 
resonance  frequency  and  quality  factor  through  calculation  of 
the  pressure  amplitude  at  the  opposite  end  as  a  function  of 
frequency.  The  result  of  the  computer  implementation  are 
compared  with  measurements  made  with  open  and  closed  ended 
prime  mover  (Ref.  2,3],  and  with  the  results  of  a  standing 
wave  analysis  of  prime  movers  [Ref.  4]. 
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IZ .  TBSORY 


The  goal  of  this  chapter  is  to  describe  the  basic 
geometry  of  a  prime  mover  and  to  summarize  the  important 
points  of  Arnott's  method.  Parallel  with  this  we  explain  how 
this  theory  is  implemented  in  the  computer  program.  The  reader 
is  referred  to  (Ref.  1]  for  full  details  of  Arnott's  method. 

A.  ARMOTT'S  METHOD. 

The  necessary  background  for  this  thesis  can  be  explained 
with  the  help  of  Fig.l  which  shows  the  basic  construction  of 
the  thermoacoustic  prime  mover.  It  constists  of  a  circular 
cross  section  acoustic  resonator  ,  rigidly  capped  at  both 
ends.  Situated  within  the  resonator  are  the  thermoacoustic 
elements  consisting  of  an  ambient  heat  exchanger,  a  prime 
mover  stack  and  a  hot  heat  exchanger.  In  our  prime  mover  both 
the  stack  and  heat  exchangers  are  parallel  plates  spaced  by  a 
few  thermal  penetration  depths.  Complete  details  of  the  prime 
mover  construction  are  given  in  [Ref.  2,3,4].  The  traverse 
coordinates  of  the  prime  mover  are  taken  to  be  x  and  y,  and 
the  longitudinal  coordinate  is  z.  The  hot  and  cold  sections 
are  circular  ducts  open  at  the  heat  exchanger  end  and  open  or 
closed  at  the  other  end.  The  prime  mover  is  divided  into 
either  5  or  7  sections  depending  on  the  boundary  conditions 
[Fig.  1,2].  The  stack  is  divided  into  11  subsections. 
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rigura  1.  Rigid  end  Prime  Mover 
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rigur*  2.  Open  end  Prime  Mover 
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The  computer  program  finds  the  normalized  acoustic 
pressure  and  the  specific  acoustic  impedance  in  all  parts  of 
the  prime  mover  as  a  function  of  frequency.  The  pressure  as  a 
function  of  frequency  gives  the  frequency  response.  From  the 
frequency  response  we  can  determine  the  quality  factor  (Q)  and 
the  resonance  frequency  (fo)  .  To  do  that  we  assume  a  linear 
temperature  gradient  across  the  stack  v’a  the  hot  and  ambient 
(cold)  heat  exchangers.  The  stack  i  divided  into  eleven 
subsections.  A  constant  temperature  is  assumed  in  each 
subsection  of  the  stack.  In  all  others  parts  of  the  prime 
mover  we  assume  that  the  temperature  is  constant  and  equal  to 
that  of  the  nearest  heat  exchanger,  ambient  temperature 
(To)  is  taken  to  be  a  function  of  z  (Tc{2))  only  in  the  stack 
area. 

The  program  starts  with  specific  acoustic  impedance  at 
the  right  end.  There  are  two  boundary  conditions  of  interest: 
rigid  and  open.  The  specific  acoustic  impedance  of  a  rigid  end 
is : 


P 


P  qC 
(1)  n 


i 


(1) 

i 


as  described  in  (Ref.  5;p.529).  In  the  above  equation  Np,  is 
the  Prandtl  number,  Cp is  the  isobaric  heat  capacity  per  unit 
mass,  Y  is  the  ratio  Cp/c,  where  c,  is  the  constant  volume  heat 
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capacity  per  unit  mass,  t)  is  the  viscosity  and  Po  represents 
the  ambient  value  of  the  density.  The  specific  acoustic 
impedance  of  an  open  end  is: 


Z*cp  oka{i^-0.6i) , 


(2) 


where  a  is  the  radius  of  the  tube  and  lc=<a/c  is  the  wavenumber, 
0)  is  the  angular  frequency  and  i=(-l)'*  [Ref.  6:p.202]  . 

The  acoustic  pressure  is  specified  at  the  same  end. 
Because  the  Q  is  independent  of  the  absolute  value  of  the 
acoustic  pressure  in  linear  theory,  the  specified  pressure  is 
only  a  relative  or  normalized  pressure.  Next  we  find  the 
pressure  and  the  impedance  at  the  entrance  to  the  hot  heat 
exchanger  using  the  Rayleigh' s  impedance  translation 
theorem  [Ref -IfS] .  This  is  one  difference  between  Arnott's  and 
other  methods.  That  is,  other  methods  find  the  fluid  variables 
by  integrating  the  governing  equations  along  the  prime  mover. 
The  impedance  translation  theorem  state  that:  within  any 
homogeneous  layer,  with  intrinsic  (or  characteristic}  specific 
acoustic  impedance  Zi„t,  the  local  specific  impedance  Z(z-l)  at 
2-1  ij  related  to  that  at  z  by  [Ref.  5:p.l39J: 


ziz-l)nz 

^"*Zj„tC08(JcJ> -lz(y)ain(kl)  ‘ 


(3) 
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Likewise,  the  acoustic  pressure  at  the  entrance  to  the  hot 
heat  exchanger  is  found  through  the  pressure  translation 
theorem: 


Pj(z-I)  «P^(z){cos  (lei)  -1  [.^1!^]  sin  (lei)),  <4) 

Z\Z) 

where  Z(z)  is  the  specific  impedance  at  2,  Pi(2)  is  the  first 
order  value  of  pressure,  i.e,  the  acoustic  pressure,  and  Zi„t 
is: 


Z.  (5) 

[QF(X  )le] 

In  the  above  equation,  Q  represents  the  porosity  and  is  equal 
to  NA,  where  A  is  the  cross  sectional  area  of  single  pore  and 
N  is  the  number  of  pores  per  unit  area,  k  is  the  complex  wave 
number  in  the  pore  and  given  by: 

(lc(X  ,X  ,))»-^.-|-^(Y  -(Y  -DPU  ,))  .  (6) 

C*  ^VA  )  ^ 


The  F(X)  and  F(Xt)  factors  arise  from  the  eq'.iation  for 
the  z  component  of  the  acoustic  velocity: 


u  ,(x,y,z) 


F(x.y;X  )  dPx<g) 
iw  p  0  dz 


(7) 
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The  F(x,y;X)  satisfies  the  following  differential  equation 
(Ref.  1): 


F(x,y;A  )  V?/*{x,y/A  )  «1,  (8) 

lA  * 

where  the  Laplacian  operator  V^j  is  given  by: 


’^^x*  0y» 


(9) 


In  our  case  where  the  pores  of  the  stack  are  parallel  plates, 
the  F  depends  only  on  y  and  the  shear  wave  number  X  or  the 
dimensionless  thermal  disturbance  number  X,  which  are  given 
by: 


and 


n 


<10) 


A 


(11) 


The  F(y;X)  is  given  by: 


F{y,A  )-l- - 1-5-.  (12) 

C08h(/^y  ) 
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In  the  above  equation  a  is  the  half  separation  distance 
between  the  parallel  plates  of  the  stack  or  heat  exchanger. 
The  pore  average  of  F(y;X)  for  heat  exchangers  and  stack,  is 
given  by: 

F(A  )  tanh{/^^)  .  (13) 

Finally  we  use  the  boundary  layer  approximation  in  the 
circular  regions  of  the  prime  mover: 

>»l-(5 /«)  (l^l) .  (14) 

In  the  above  relation  R  is  twice  the  hydraulic  radius  of  the 
pore,  or  twice  the  ratio  of  the  transverse  pore  area  to  the 
pore  perimeter.  Actually  for  heat  exchangers  and  stack  with 
parallel  plates,  R  is  the  separation  distance  (2a)  between  the 
plates  and  for  the  circular  parts  it  is  the  radius  of  the 
tube.  $  is  the  viscous  or  thermal  penetration  depth. 

Once  Pi  and  Z  are  known  at  the  entrance  to  heat  exchanger, 
we  need  the  values  just  inside  the  heat  exchanger.  Using  the 
pressure  and  the  impedance  translation  theorem  we  have  a 
solution  up  to  the  stack  where  a  linear  temperature  gradients 
exists.  In  the  heat  exchanger  three  things  are  different  from 
the  circular  parts.  These  are:  the  wavenumbers  in  the  heat 
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exchangers  are  complex,  the  porosity  now  is  NA,  and  F(X)  as  we 
mention  above  is  given  by  equation  (13)  . 

The  point  now  is  how  to  find  pressure  and  impedance  in 
the  stack.  In  this  case  Arnott's  method  uses  the  following 
system  of  the  first  order  differential  equations  that  can  be 
solved  using  Range-Kutta  methods: 

■^2M-iJc(z)Z^e(z)  (1 — (z)Zlz),  (15) 
dZ  Zlattz)* 

and 


dP^{2) 

dz 


•ik(z)  Zj,nt(z) 


Pi(z) 

Z(z)  ' 


(16) 


where  a(z)  is: 


a(A  .X  ^ 


ZILeL-t 

(F(A  ))  \ 
1-N^ 


(17) 


To,  is  the  temperature  gradient  given  by: 


r  ,dT„(z) 
dz  • 


(18) 


P  is  the  thermal  expansion  coefficient: 
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(19) 


.(ijBL 
‘  e  T 


^  p 


In  our  program,  this  system  of  equations  is  solved  with  the 
Runge-Kutta-Fehlberg  method  for  every  subsection  of  the  stack. 
In  so  doing,  the  impedance  and  the  pressure  distribution  in 
the  stack  i.s  obtained.  Now  using  the  pressure  and  the 
impedance  translation  theorem  we  find  the  pressure  and  the 
impedance  in  the  rest  of  the  prime  mover. 

Knowing  the  pressure  at  the  left  end  we  can  calculate  the 
quality  factor,  the  resonance  frequency  and  the  maximum 
pressure  amplitude.  This  is  done  by  minimizing  the  following 
equation  (Ref.  8]: 


.SL 


f  f. 


0* 


■•‘aaplic*. 


(20) 


where  Q  is  the  quality  factor,  |P(1)|  is  the  maximum 
amplitude  of  pressure  at  the  left  end,  fo  is  the  resonance 
frequency,  amplit  is  the  calculated  amplitude  of  pressure  at 
the  left  end  and  is  a  function  of  the  frequency  f.  In  the 
above  relation  there  are  two  known  variables  f  and  amplit  and 
three  unknowns  (Q,  fo/P(l))  that  we  want  to  find.  To  find  these 
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unknowns,  the  3-parameter  least  squares  technique  is  used  with 
the  help  of  simplex  method. 

Up  to  this  point,  Q  and  have  been  determined  only  by 
considering  the  pressure  amplitude  at  the  left  end.  We  have 
not  worried  about  matching  the  acoustic  impedance  at  the  left 
end.  Impedance  matching  allows  us  to  refine  the  values  of  Q 
and  fa  through  the  complex  eigenfrequency .  Using  this 
technique,  one  assumes  the  angular  frequency  is  complex  and 
given  by: 


<21) 

Substituting  the  initial  value  of  f^  and  Q  determined  from  the 
frequency  response,  into  this  equation  gives  an  initial  guess 
at  the  complex  eigenfrequency  for  the  system.  This  frequency 
is  used  to  calculate  the  impedance  throughout  the  prime  mover 
as  explained  earlier.  The  value  of  impedance  in  the  gas  at  the 
left  end  is  compared  to  the  value  for  the  boundary  condition 
(i.e.  the  impedance  of  a  rigid  end  with  thermal  loss).  The 
difference  between  the  computed  and  actual  impedance  is  used 
to  adjust  the  eigenfrequency.  The  impedance  is  computed  with 
the  new  eigenfrequency  until  both  real  and  imaginary  parts  of 
the  impedance  matche  at  the  left  end.  The  final  values  of  fo 
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and  Q  are  determined  from  the  real  and  imaginary  parts  of  the 
eigenf requency  that  best  matche  the  impedance. 


B.  IMPLEMENTATION  OF  ABNOTT' S  METHOD 

The  program  uses  functions  that  already  exist  in  Matlab 
where  possible.  For  the  Runge-Kutta-Fehlberg  5th  order  method 
we  modified  the  existing  function  in  Matlab  to  integrate  with 
initial  values  of  the  integration  variable  greater  than  the 
final.  This  is  necessary  because  we  put  the  origin  at  the  left 
end  (2=0)  and  we  begin  our  calculations  from  the  right  end 
(z=L) .  One  reason  that  we  use  this  particular  Runge-Kutta 
method  is  its  accuracy.  The  Runge-Kutta-Fehlberg  computes  two 
Runge-Kutta  estimates  for  the  value  y„,i  but  of  different 
orders  of  error.  The  global  error  in  this  numerical  method  is 
O(h')  and  the  local  error  is  0(h‘)  iRef.  7). 

To  do  the  minimization,  we  use  the  function  fmins  with 
the  help  of  simplex  numerical  method.  The  fmin  function  with 
simplex  method  is  used  to  match  the  impedances  at  the  left  end 
using  the  technique  of  least  squares. 
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ZII.  THE  METLXB  PROGRAM 

The  goal  of  this  chapter  is  to  provide  an  overview  of  the 
structure  of  the  Mat lab  program  that  implements  Arnott's 
method.  Furthermore  it  describes  all  the  functions,  special 
tricks  and  limitations  of  this  code.  It  also  explains  how  to 
run  the  program. 

A.  OOTLIHS  or  THE  MATLAB  PROGRAM 

The  program  is  divided  into  three  parts.  In  the  first 
part  the  variables  are  specified  for  every  section  of  the 
prime  mover.  They  are:  the  type  of  section  (open,  heat 
exchanger,  stack) ,  the  number  of  subsections  that  exists 
inside  the  section,  the  temperature  at  the  left  and  right  ends 
of  the  section,  twice  the  hydraulic  radius  of  the  section,  and 
the  porosity  of  the  section.  Also,  specified  are  the  frequency 
range  and  the  number  of  intervals  in  this  range. 

The  second  part  is  the  main  program.  It  finds  the 
impedance  and  the  pressure  distribution  in  every  section  of 
the  prime  mover  for  each  frequency  within  the  assumed  range. 
To  accomplish  this  it  begins  at  the  right  end  of  the  tube 
where  it  calculates  the  impedance  of  the  end,  rigid  or  open, 
and  assigns  an  arbitrary  value  to  the  pressure.  Continuing,  it 
calculates  using  the  translation  theorems,  the  pressure  and 
the  impedance  in  all  other  sections  except  in  the  stack.  In 


the  stack,  where  there  exists  a  temperature  gradient,  the 
program  calls  the  function  "trode45”  (a  modified  version  of 
the  function  of  Mat  lab  *’ode45")  to  integrate  from  the  right 
end  to  the  left  end  of  the  stack.  This  function  solves 
differential  equations  using  the  5th  order  Runge-Kutta 
numerical  method  with  Fehlberg  coefficients. 

The  last  part  of  the  Matlab  program  uses  the  calculated 
pressure  at  the  left  end  from  the  previous  part  to  find  the 
maximum  amplitude,  the  resonance  frequency  and  the  quality 
factor.  These  parameters  are  found  from  a  least  square  fit  to 
equation  (20)  .  The  minimization  is  performed  with  the  matlab 
function  "fmins"  with  "options  1,2,3,14".  Finally,  the  Q  and 
fo  estimates  are  used  as  a  starting  points  for  the  final 
calculation  of  the  quality  factor  and  the  resonance  frequency 
using  the  complex  eigenfrequency  method.  This  refinement 
usually  is  a  small  correction  (in  the  fourth  decimal  place)  to 
initial  guesses.  Finally,  the  program  plots  1/Q  and  f^  versus 
temperature  difference  and  saves  the  results  in  a  data  file. 


B.  LZNXTATKHIS  OF  TBB  FROGRMf,  SPKCZAL  POZMTS. 

There  are  two  programs,  one  for  the  open  end  prime  mover 
and  one  for  the  rigid  end  tube.  To  run  the  open  end  progam, 
one  needs  to  run  "arffl.m"  with  the  functions  "error,m- 
erro3.m-argo3f .m" .  To  run  the  rigid  end  program  one  needs  to 
run  "argo2.m"  with  the  functions  "error. m-erro3. m-argo3. m" . 
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The  difference  in  the  two  progams  is  that  one  uses  "argo3.m" 
and  the  other  uses  "argof3.m”.  Before  running  one  program  or 
the  other,  one  needs  to  switch  these  functions  in  the 
"erro3.m"  function. 

Futhermore,  there  are  some  factors  that  influence  the 
operation  of  the  program.  One  important  step  is  the  choice  of 
the  initial  range  of  the  frequencies.  To  begin  the  program  one 
needs  to  know  the  approximate  resonance  frequency  as  well  as 
which  longitudinal  acoustic  mode  is  to  be  used.  Initially  the 
frequency  range  is  fo  ±  20Hz.  After  the  program  runs  for  the 
first  value  of  AT  the  frequency  range  is  automatically 
adjusted  to  be: 


Another  important  factor  is  of  course  the  number  of 
points  between  the  maximum  and  minimum  frequency  that  the 
program  uses.  Usually  this  number  is  500  but  sometimes  when 
the  values  of  quality  factor  gets  large  this  may  have  to  be 
doubled  (1000  points) .  In  these  cases,  it  may  also  be 
necessary  to  use  steps  of  5  degrees  Kelvin  instead  of  the 
usual  10  degrees  Kelvin. 

Finally,  another  important  point  concerns  the  "counters" 
that  are  used  by  the  program.  The  most  important  counter  is 
the  "flag".  If  somebody  wants  to  run  the  program  for  a  case 
where  the  quality  factor  increases  with  AT,  one  needs  to  put 
the  initial  value  of  the  "flag"  counter  equal  to  zero. 
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otherwise,  it  should  be  set  equal  to  one.  The  purpose  of  this 
counter  is  to  force  the  program  to  scan  for  negative  values  of 
the  inverse  quality  factor  at  the  proper  time. 
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XV.  VALIDATION  AND  RESULTS 

In  this  chapter  we  compare  the  results  of  the  program  to 
experimental  results  for  two  different  prime  movers  [Ref.  2 
and  3]  as  well  as  with  the  results  of  the  standing  wave 
analysis  by  Atchley  [Ref  4]  . 

The  first  results  are  for  a  rigid  end  prime  mover 
described  in  Ref.  2  and  4.  Comparison  of  both  1/Q  and 
resonance  frequency  are  made.  Figures  3  through  5  show  1/Q  vs 
AT  (in  degrees  K)  for  three  different  mean  gas  pressures.  The 
gas  is  helium.  In  these  figures,  the  open  circles  represent 
measured  values,  the  *'s  represent  the  results  of  the  standing 
wave  analysis  and  the  solid  line  represents  the  result  of  the 
Matlab  program.  The  overall  agreement  between  the  present 
analysis  and  the  measurements  is  good,  but  not  as  good  as  for 
the  standing  wave  analysis.  The  present  analysis  agrees  best 
at  low  values  of  AT.  The  reasons  for  the  increased  discrepancy 
at  higher  values  of  AT  are  not  understood.  One  would  expect 
the  present  analysis  to  agree  better  with  measurements  than 
the  standing  wave  analysis.  This  is  expected  because  the 
present  analysis  makes  fewer  assumptions.  In  the  particular, 
it  does  a  better  job  of  tak'.ng  into  account  changes  in  the 
acoustic  pressure  amplitude  in  differents  parts  of  the  prime 
mover.  The  standing  wave  analysis  ignores  any  changes  in 
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pressure  amplitude  due  to  changes  in  cross  sectional  area  of 
the  prime  mover. 

Figures  6  through  8  show  results  of  resonance  frequency 
(in  H2)  vs  at  for  the  same  three  mean  pressures.  The  agreement 
between  the  present  analysis  and  the  measured  frequencies  is 
within  1%  in  all  cases.  The  standing  wave  analysis  also  agrees 
with  the  measured  values  to  about  the  same  accuracy.  However, 
the  two  methods  show  different  depedences  on  AT.  This  is 
because  1)  the  standing  wave  analysis  considers  only  the 
temperature  dependence  of  the  sound  speed,  whereas  the  present 
analysis  includes  the  effects  of  dispersion  due  to  the 
boundary  layers,  and  2)  the  standing  wave  analysis  ignores  the 
finite  impedance  of  the  rigid  ends.  Of  these  two  differences 
the  first  is  the  more  important.  This  results  in  different 
curvatures  of  the  two  theoretical  lines. 

The  second  set  of  comparisons  is  made  for  an  open  ended 
prime  mover  used  in  Che's  thesis  research  (Ref.  3].  The 
unusual  property  of  this  prime  mover  is  that  while  1/Q 
decreases  with  AT  for  the  first  and  third  modes  (as  is  usual 
for  a  prime  mover),  1/Q  increases  with  AT  for  the  second  mode. 
This  behavior  provides  a  good  test  for  theoretical  models.  The 
results  for  1/Q  vs  AT  are  shown  in  figures  9  through  11  for 
the  first,  second  and  third  modes,  respectively.  It  is  seen 
that  the  present  analysis  and  the  standing  wave  analysis  are 
in  close  agreement  in  the  first  and  second  modes.  The 
agreement  diminishes  slightly  at  the  third  mode.  The  agreement 
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with  the  measured  values  is  worse  than  before.  However,  this 
is  most  likely  due  to  contamination  of  the  helium  with  air  and 
water  vapor  as  discuseed  in  Che's  thesis.  Another  interesting 
aspect  of  the  analysis  is  the  prediction  of  an  extremum  in  1/Q 
for  the  second  and  third  modes. 


rigur*  3.  Graph  of  1/Q  vs  AT  (K)  for  the  closed  end  prime 
mover  filled  with  helium  at  170  kPa  mean  pressure. 


21 


FrequetM:y 


Negative  DeltaT 

Figure  6.  Graph  of  f„  (Hz)  vs  AT  (K)  for  the  closed  end  prime 
mover  filled  with  helium  at  170  )cPa  mean  pressure. 
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mover  filled  with  helium  at  500  hPa  mean  pressure. 
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rigur*  9.  Graph  1/Q  vs  AT  (K)  for  the  first  mode  of  the  open 
ended  prime  mover  filled  v'ith  helium  gas  at  101  kPa. 


open  ended  prime  mover  filled  with  helium  gas  at  101  kPa. 


/ 

l 


rigura  13.  Graph  of  fo  (Hz)  vs  AT  (K)  for  second  mode  of  the 
open  ended  prime  mover  filled  with  helium  gas  at  101  )cPa. 


Resonance 


the  open  ended  prime  mover  filled  with  helium  gas 
at  101  kPa. 
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V.  soioaRy 


The  results  of  this  thesis  can  be  summarized  as  follows. 
Arnott's  method  was  implemented  in  a  Matlab  progam.  The 
impedance  and  the  pressure  distribution  were  calculated  along 
the  prime  mover.  From  the  frequency  response  at  the  left  end 
and  complex  eigenfrequency,  the  quality  factor  and  the 
resonance  frequency  was  estimated  for  both  cases  rigid  and 
open  end  prime  movers.  The  results  of  the  analysis  were 
compared  with  those  of  a  standing  wave  analysis  and  with 
experimental  results.  In  the  most  of  the  cases  the  results 
agree  well.  A  future  investigation  can  be  the  examination  of 
the  output  work  of  the  prime  mover. 
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APPIMDZX  A.  THE  NATLAB  PROdUM  POR  TBB  OPKM  BHD  PRIMS  MOVER 


%  Program  for  prime  mover. 

%  PART  1  (CHANGE) 

Clear; 

clg; 

coun=0;  %  counter 

pt=0;  %  counter 

qual=2eros (45, 1) ;  %  quality  factor 

inqual=2eros (45, 1) ;  %  inverse  quality  factor 

frO ’2eros (45, 1) ;  %  resonant  frequency 

deltaT^seros (45, 1) ;  %  temperature  difference 

point =0;  %  indicator 

flag=l;  %  indicator 

slope=sum (' negative' ) ; 

load  open.dat  %  standing  wave  analysis  data 

%  Main  Loop  for  evrey  temperature  difference  (itrgh) 
for  itrgh»293;-*5: 158 
coun*l+coun 

numfre>1000;  %  number  of  frequencies 

%  establish  some  often  used  constants 

s«5  ;  %  number  of  sections 

j=sqrt (-1) ; 
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npr=2/3* (1+eps) ; 

ganiina*5/3*  (1+eps)  ; 

cp* (2. 5*8. 3143/ (4.e-3) )* (1+eps) ; 

ksolid=. 16* (1+eps) ; 

sqri= ( (1+ j) /sqrt (2) ) * (1+eps) ; 

pl*ones(lb,numfre) .*(10*(-8));  %  presure  inside  the  tube 
f req^zeros (1, numf re) ;  %  frequncies 

w=zeros (1< numfre) ;  %  angular  frequencies 

teren=zeros (1/ s) ;  %  end  of  sections 

numblay=zeros (l,s) ;  %  layers  of  sections 

lengsec»zeros (1, s) ;  %  lengths  of  sections 

tempR*zeros (1, s) ;  %  rigth  temperature  of  sections 

tempL»zeros (1, s) ;  %  left  temperature  of  sections 

poros*zeros  (1,  s) ;  %  porocity  of  stac)c 

ratio«zeros (1, s) ;  %  ratio 

%  SET  VALUES 

fremin-864;  %  minimum  frequency 

fremax»904;  %  maximum  frequency 

termin=sum (abs (' free' ) ) ;  %  end  of  the  last  section 

ampres-(1.01e+5)*  (1+eps);  %  ambient  presure 

dramp*(l.e-l6);  %  driver  pressure 

%  section  "1" 

terenl=' opentu' ; 

teren (1,1) »sum (abs (terenl) ) ; 

numblay (1, 1)»1; 
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lengsec (1, 1) * (66.25e-2) * il*eps) 
tempRd,  1)  >293*  (1-feps)  ; 
teinpL(l,l)*293*  (1+eps); 
ratio  (1, 1)  *  (1 . 917e-2)  *  (1-t-eps)  ; 
poros (1, 1) =1* (1+eps) ; 

es(l,  :)»pi*ratio(l, :)  *  (1+eps) 
%  section  "2" 

teren2=' hexch' ; 

teren (1, 2) >sum (abs (teren2) ) ; 

numblay (1, 2) >1; 

lengsec (1, 2) > ( . 82e-2) * (1+eps) ; 
terapR(l,2)=293* (1+eps); 
tempi (1, 2) >293* (1+eps) ; 
ratio  (1, 2) ■ <5.08e-4) * (1+eps) ; 
poros (1, 2) ». 666* (1+eps) ; 

%  section  "3" 
teren3*'  stac)c' ; 
teren (1, 3) »sum(abs (teren3) ) ; 
numblay (1, 3) >11; 
lengsec (1, 3)« (2.59e-2) * (1+eps) ; 
tempRd,  3) -itrgh*  (1+eps) ; 
tempi (1, 3) >293* (l*eps) ; 
ratio (1, 3) ■ (8. 636e-4) * (l+eps) ; 
poros (1, 3) 89473* (1+eps) ; 

%  section  "4" 


teren4-*hexch' ; 


teren (1, 4) *sum (abs (teren4) ) ; 
numblay (1, 4) =1; 

lengsec (1, 4) » ( .82e-2) * (1+eps) ; 
tempRd,  4)  =itrgh*  (1+eps)  ; 
tempLd,  4)  »itrgh*  d+eps)  ; 
ratiod,  4)  =  (5.086-4)  *  (l+eps); 
poros  (1,  4)  *. 666*  d+eps) ; 

%  section  "5" 

teren5=' opentu' ; 

teren (1, 5) -sum(abs (terenS) ) ; 

numblay (1, 5) »1; 

lengsec (1, 5) * ( . 6852) *  (1+eps) ; 

tempRd,  5)  »itrgh*  (1+eps)  ; 

tempL(l,5)=itrgh* (1+eps); 

ratio  d,  5) » (1 . 917e-2) * (l+eps) / 

poros  (1, 5)  *1*  d+eps) ; 

%  PART  2 
%  MAIN  PROGRAM 
%  FREQUENCY  DISTRIBUTION 

numtot^suin  (numblay)  +1; 

pt*pt+l; 

if  pt-*l 

df*fr0  (coun-* ) /abs (qual (coun-1) ) 
f remin-frO (coun-1) -df; 
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f remax* frO  <coun-l) ♦df ; 

end 

i*l :numfre; 
if  i«-l 

f req*f remin; 
else 

freq  (1/  i)  *f remin »■  (fremax-f  remin)  .  *  (i .  /  (numfre-l) ) ; 

end 

w*2*pi*freq; 

%  GET  THE  SPECIFIC  ACOUSTIC  IMPENDANCE  AND  PRESSURE  AT  ALL 
POINTS. 

%  START  AT  THE  RIGHT  AN  .MOVE  AT  THE  LEFT, 
dens  (l,numtot)  *ampres*4.0e-3  /  (tempRd,  s)  *8. 3143) ;  %  density 
vise  (1,  numtot)  *1 . 887e-5*  (tempRd,  s)  /273. 15)  *  ( .  6567) ;% 

viscosity 

kgas (1, numtot) “Vise d, numtot ) *cp/npr;  %  termal  conductivity 
sspeedd, numtot)  *972. 8*sqrt  (tempRd,  s) /273. 15) ;  %  speed 

%  isothermal 

typel*' free' ; 
type2*' rigid' ; 

if  termin**sum<abs (type2) ) 

%  IMPEDANCE  OF  RIGID  TERMINATION. 

fac (numtot, 1 ;numfre) «sqrt ( (dens (1, numtot) . . . 

* s speed (1, numtot) *2) . / (w. * vise (1, numtot) ) ) ; 
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z  (nujiltot/ 1  rnumfre)  ■  (1+ j)  •  *  (dens  (1#  numtot)  .  *sspeed  (1,  numtot )  . 
. . .  *fac  (numtot,  1  :nuinfre)  *sqrt  (npr) )  /  (sqrt  (2)  *  (gaimna-l) ) ; 
%  IMPEDANCE  OF  FREE  TERMINATION 
else 

z  (numtot,  1  :nunif  re)  »free  (dens  (1,  numtot) ,  vise  (1,  numtot) ,  w,  rati 
0(1, s) , s speed (1, numtot) , gamma, npr) ; 
end 

%  IMPEDANCE  TRANSLATE  FOR  THE  OPEN  TUBE  OR  HEAT  EXCHANGER, 
fault*'  stac)c' ; 
truel*' opentu' ; 
true2*'hexch' ; 
for  )c=s:“l:l 
jup«0; 

jup»numtot-suro(  jup+numblay  (l,)c:s) )  ; 
if  teren  (1,  )c)  *»sum(abs  (fault) ) 
dens  (1,  jup)  =ampres*4 .Oe-3  /  (tempR(l,)c) .  *8.3143) ;  %  density 
vised,  jup)»1.887e-5.*  (teropR(l,)c)  /273.15)  .*  (.6567)  ; 

%  viseosity 

)c9as  (1,  jup)  *vise  (1,  jup) .  *ep/npr;  %  termal  eonduetivity 
sspeedd,  jup)  *972. 8.  *sqrt  (tempRd, k) /273. 15)  ;  %  speed 

%  isothermal 

%  GET  LAMBDA  AND  LAMBDA (T) 

lambda (jup, 1 :numfre) *ratio (1, k) . *sqrt (dens (1, jup) . *w. /vise (1 
»  jup) ) ; 

lambdt (jup, 1 :numfre) -sqrt (npr) . *lambda ( jup, 1 :numfre) ; 


%  GET  F(L),  F(L(T))  AND  WAVENUMBERS  FOR  OPEN  TUBE 

PARTS 

if  teren  (1,  k)  »*suni<abs  (truel) ) 
flam(jup,  l:numfre)«l-(l+j)  .*sqrt  (2)  . /lambda  (jup,  1  :nuinf  re) ; 

%  f(l) 

flamt ( jup, l:numfre)»l-(l+j) .*sqrt (2) . /lambdt ( jup, linumfre) ; 

%  f(l(t)) 

facl» (1+ (gamma-l) /sqrt (npr) ) /sqrt (2) ; 
ka ( jup, 1 :numfre) » (w./sspeed(l, jup) ) . . . 

(1+  (1+j) .*facl. /lambda (jup, l:numfre) ) ; 
else 

%  GET  F(L),  F(L(T))  AND  WAVENUMBERS  FOR  HEAT 
EXCHANGERS  TYPE  "SLIT". 

sqrmi» (1- j) /sqrt (2) ; 

ar ( jup, Irnurof re) -sqrmi. "lambda (jup,  l:numfre) ;  / 

ar9um(jup, l:numfre)>exp(-ar(jup, l:numfre) ) ; 
ar ( jup, 1 : numf re) «ar ( jup, 1 : numf re) /2 ; 
ctanh( jup, 1 :numfre) « (l-argum( jup, l:numfre) ) . / (l+argum(  jup, 1: 
numf re) ) ; 

flam  (jup, 1 : numf re) *l-ctanh (jup, l:numfre) ./ar ( jup, Itnumfre) ; 

%  f(l) 

ar (jup, 1: numf re) "sqrmi* lambdt (jup, 1: numf re) ; 
argum( jup, 1 : numf re) »exp(-ar ( jup,  1 : numf re) ) ; 
ar ( jup, l:numfre)*ar ( jup, 1: numf re) /2; 
ctanh(  jup,  1  rnumfre) * (l-argum(  jup,  l;numfre) ) ./  (l-t'argum(  jup,  1; 
numf re) ) ; 
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flamt ( jup/ 1 :numfre) =l-ctanh ( jup,  1 :numfre) . /ar ( jup, 1 :numfre) ; 

%  f{l{t)) 


low=w/sspeed{l, jup) ; 

ka (jup, 1 cnumfre) =low. *sqrt ( (gamma- (gamma-1) . . . 

. lamt (jup, Itnumfrc) ) . /flam( jup, 1 :numfre) ) ; 
end 

%  TRANSLATION  THEOREM 

comden ( jup, 1 rnumfre) «dens (1,  jup) . /flam (jup, 1 :numfre) ; 
zint (jup, 1 :numf re) =comden ( jup,  1 :numf re) . *w. / (ka ( jup, 1 :numfre 
) . *poros (1, k) ) ; 

dsub (1, k) »lengsec (1, k) . /numblay (1, k) ; 
ar( jup, 1 :numfre) -ka ( jup,  l:numfre) .*dsub(l,k) ; 
sn( jup, l:numfre)»sin(ar( jup, l:numfre) ) ; 
cs ( jup, l:numfre)»cos (ar ( jup, 1 :numfre) ) ; 
ct  (jup,  l:numfre)«cs  (jup,  l:numfre)  ./sn(  jup,  l:nuinfre)  ; 
fac3  ( jup,  1  :numf  re)  a:zint  (jup,  1  :numfre)  .  /z  ( jup-fl,  1  :numfre) ; 

z  (jup,  1  :numfre)  «zint  (jup,  1  rnumfre) .  *  (ct  (jup,  1  rnuaifre) . . . 
- j . *fac3 (jup, 1 :numfre) ) . / (fac3 ( jup, 1 :numfre) . *ct (jup, 1 :numfr 
e)-j)  ; 

pi (jup, 1 znumfre) »pl ( jup+1, 1 rnumfre) . * (cs ( jup, 1 :numfre) . . . 
- j . *fac3 ( jup, 1 rnumfre) . *sn ( jup, 1 inumfre) ) ; 
else 

%  STACK  WITH  LINEAR  DISTRIBUTION  OF  TEMPERATURE  BETWEEN 
THE  ENDS  OF  THE  STACK 

toz=(tempR(l,k)-teropL(l,k) )/lengsec(l,k) ;  %  approximation 

dz/dt 
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ns»0; 

for  la*k+numblay(l,k)-l:-l:)c 
ns»ns+l; 

tens*tempR(l,)c)-(tempR(l,)c)-tempL(l,k) )  *ns/nujnblay  (1,  k) ; 
tnsml=tempR(l,  k)  -  (tempRd,  k)  -tempL  (1,  k) )  *  (ns-1)  /numblay  (1,  k) 

t 

tave (1, la) * (tens+tnsml) /2; 
dsubd,  k) slengsec  d, k)  /numblay  (1, k) ; 

densd,  la)»ampres*4.0e-3  /  (tave  (1,  la)  *8 . 3143) ;  %  density 

vise (1, la) =1 . 887e-5* (tave d, la) /273. 15) * ( . 6567) ;  %  viscosity 
kgas (1, la) =visc (1, jup) *cp/npr;  %  termal  conductivity 
sspeedd,  la)«972.8*sqrt  (taved,  la) /273.15) ;  %  speed 

%  isothermal 

%  GET  LAMBDA  AND  LAMBDA (T) 

lambda  (la,  1  :numfre)  ^ratio  (a,  k) .  *sqrt  (dens  (1,  la)  .  '*w.  /vise  (1, 1 
a) )  ; 

lambdt (la, 1 :numf re) «sqrt (npr) ^lambda (la, 1 :numfre) ; 
ar (la, 1 :numfre) ^sqrmi* lambda (la, Itnumfre) ; 
argumda,  1  :numfre)  sexp(-ar  (la,  1  :numfre) ) ; 
ar (la, 1 :numfre) ^ar (la, l:numf re) /2; 
ctanhda,  1  znumfre) > (l-argumda,  1  :numfre) ) . /  (1  ■•■argumda,  1  :num 
fre) ) ; 

flam (la, 1 :numfre) ^l-ctanh (la, l:numfre) . /ar da, 1 :numfre) ; 

%  fd) 

ar (la, 1 :numfre) =sqrmi*lambdt (la, 1 inumfre) ; 
argumda,  1  :numfre)  "exp (-ar  (la,  l:numfre) ) ; 
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ar (la, 1 :nurof re) =ar (la, 1 :numf re) /2; 

ctanhda,  1  .-numfre)  =  (l-argumda,  l:nuinfre) )  ./  (l+argum (la,  l:num 
fre) ) ; 

f lamt  da, 1 tnumfre) =l-ctanh (la, 1 :numfre) . /ar (la, 1 inumfre) ; 

%  f(l(t)) 

lowsw/sspeecKl,  la) ; 

)ca  (la,  1  :numf  re)  =low.  *sqrt  ( (gaituna-  (gamma-1)  . . . 

. *flamt (la, 1 :numfre) ) . /flam (la, 1 :numfre) ) ; 
zint (la, 1 :numf re) *dens (1, la) . *w. . . 

.  /  (poros  (1,  )c)  .  *  flam  (la,  1  :numf  re)  .  *)ca  (la,  1  :numfre) ) ; 
alprimda,  1  :numfre)  »toz*  (flamt  (la,  1  :numfre) . . . 

. /flam (la, 1 cnumfre) -1) . / (2*tave (1, la) * (l-npr) ) ; 

)cal»)ca  (la,  1  :numf re) . * ; 
zintl*2int  da,  l:numfre) .  * ; 
alprl»alprim(la, 1 :numfre) . ' ; 

zetaO»  (lengsec  d,]c)  -  (ns-1)  *lengsec  (1,  Jc)  /numblay  (1,  k) ) ; 
zetafina (lengsec (1, k) - (ns) *lengsec  d, k) /numblay (1,  k) ) ; 
zO»z (la+1, 1 :numfre) . ' ; 
tol»l .e-4; 

(zeta, zout ] »trode45 (' tube' , zetaO, zetafin, zO, tol, kal, zintl, al 
prl)  ; 

n« length (zeta) ; 

z (la, 1 :numfre) »zout  (n, : ) ; 

zdumy=zO; 

plO»pl (la+1, 1 :numfre) . ' ; 
tol=l .e-4; 
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[zeta,plout  }*trodie45  (' tubel' ,  zetaO,  zetaf  in,plO/  tol,  leal,  zintl 
,  zdumy) ; 

nl>length (zeta) ; 

pi (la, 1 rnumfre) *plout (nl, : ) ; 

end 

end 

end 


%  MINIMIZATION  TO  GET  RESONANCE  FREQUENCY  AND  QUALITY 
FACTOR. 

tfun«- j*w.  *z  (1, 1  rnumfre)  *drainp.  /pi  (1, 1  rnumfre)  ; 
for  a«l:s 

pi  (a,  l;nuinfre)»pl  (a,  1  :numfre)  .*tfun; 
end 

amplit^zeros  (l,nuinfre) ; 

anpllt  (1,  l:n\iinfre)«abs  (pi  (1,  l:numfre) ) ; 

(maxamp, i ] «max (amplit ) ; 
resfre»freq(i) 
amphaf^maxarop/sqrt (2) ; 
q»0; 

frehaf»0; 
if  point"0 

for  s*2;numfre 

if  amplit (s) >«amphaf  &  amplit (s-1) <«amphaf 
frehaf= (freq(s) +freq (s-1) ) /2; 
q«0.5*resfre/ (resfre-frehaf ) 
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end 


end 

else 

for  s»2:nuinfre 

if  anplit  (s)  <>ainphaf  &  amplit  (s-1)  >«amphaf 
ama=l 

frehaf=(freq(s)+freq(s-l) ) /2 
q=0.5*resfre/ (resfre-frehaf ) 
end 

end 

end 

if  q~=0 

xOszeros (3, 1) ; 
xO  (1, 1)  «inaxamp; 
xO  (2, 1) *resfre; 
x0(3,l)»q; 
options  <1) =0; 
options (2) =1 .e-4; 
options  <3) *1 .e-6; 
options (14) >2000; 

est>fmins (' erro' ,x0, options, [ ] , amplit, freq) ; 
inqual  (coun,  1)  »est  (3, 1)  (-1) ; 
frO (coun, 1 ) =est (2, 1 ) ; 
delta! (coun, 1) >abs (itrgh-293) ; 
qual (coun, 1) >est (3, 1) ; 
end 
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if  coun-=l 


dif=inqual (coun-1, 1) -inquai (coun, 1) ; 
if  dif<0  &  coun==2 

slope=sum(abs ('positive' ) ) ; 
elseif  dif<0  &  coun=«3 

slope=sum(abs ('positive' ) ) ; 

end 

if  flag==0  &  slope~=sum(abs ('positive' ) ) 
if  dif<=0  I  q==0 
point®  L; 
f lag=l; 

for  s=2;numfre 

if  amplit (s)<=amphaf  &  amplit (s-1) >“amphaf 
frehaf® (freq(s) +freq(s-l) ) /2; 
q=0.5*resfre/ (resf re-frehaf ) 
end 
end 

xO  ( 1 , 1 )  =inaxa!np; 
xO (2, 1) =resfre; 
x0(3,l)=q; 
options (1) ®0; 
options (2) *1 .e-4; 
options (3) =1 .e~6; 
options (14) =2000; 

est=fmins ('erro' , xO, options, t ] , amplit, freq) ; 
inquai  (coun,  l)»est  (3, 1) ''  (-1) ; 
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f rO (coun, 1) =est  (2, 1) ; 
delta? (coun, 1) =abs (itrgh-293) ; 
qual (coun, 1) =est (3,1); 

end 

end 

end 

%  PART  THREE 

%  MINIMIZATION  USING  COMPLEX  ANGULAR  FREQUENCIES, 
options (1) =0; 
options (2) =1 .e-4; 
options (3) =1 .e-4; 
options  (14) *700; 

wl*2*pi*fr0 (coun, 1) * (1- j/ (2*qual (coun,  1) ) ) ; 
y2=wl+wl*10 .e-6; 
yl=wl-wl*10.e-6; 

ella=fmin (' erro3' , yl,y2, options, itrgh) ; 
frO  (coun, 1) =real (ella) / (2*pi) ; 

qual (coun, 1) »-real (ella) / (2*imag (ella) ) ; 
inqual  (coun,  1)  =qual  (coun,  1) ''  (-1)  ; 
end 

%  PLOTS  AND  STORE  THE  DATA  (CHANGE) . 

y* I (inqual (1 :coun, 1) ) ' ;  (delta? (1 :coun, 1) ) ' ; 
(frO (1 : coun, 1) ) ' ] ; 

f id=fopen (' inqf3.txt' , ' w+' )  ; 
fprintf (f id, ' Inqual=%-9. 6f  T=%-4.1f  RF=%-5.2f\n' ,y) ; 


47 


load  compf2.dat 
figure 


plot (delta? (1 :coun) , inqual (1 :coun) , compf2 ( ; ,  1) , compf2 (: , 2) , 
o'  . . . 

, open ( : , 1 ) , open ( : , 4 ) , ' * ' ) 
title (' Inverse  Quality  Factor  V.S  Temperature 
Difference. Open  End'); 

grid 

xlabel ('Negative  Delta  "T"  '); 
y label ('1/Q'); 

gtext ('Ambient  pressure  101  Kpa' ) 

gtext (' Experimental  Result  "o"') 
gtext (' Computer  Result 
gtext ('Third  Mode') 
print  -deps  inqufS 
figure 

plot (delta? (1 :coun) , frO (1 :coun) ) 
xlabel ('Negative  Delta  "T"  '); 
ylabel ('Resonance  Frequence' ) ; 
title (' Resonance  Frequency  V.S  Temperature 
Difference. Open  End'); 
gtext ('Ambient  pressure  101  Kpa') 
gtext ('Third  Mode') 
grid 

print  -deps  frf3 


%  FUNCTION  TRODE45.M 
function  . . . 

[tout  /  yout )  =trof*e45  (FunFcn,  tO,  t final,  yO,  tol/  kal,  zint  1,  alprl) 
alpha-(l/4  3/8  12/13  1  1/2]'; 

beta»[  (  1  0  0  0  00)/4 

[3  9  000  0)/32 

[  1932  -7200  7296  0  0  OJ/2197 

[  8341  -32832  29440  -845  0  OJ/4104 

[-6080  41040  -28352  9295  -5643  OJ/205201'; 

gamma«[  [902880  0  3953664  3855735  -1371249  2770201/7618050 
[  -2090  0  22528  21970  -15048  -27360} /752400 

]'; 

%  initialization 
t*t0; 

hniax>  (t-tf  inal)  /5; 
hmin* (t-tfinal) /2000; 

h» (t-tfinal) /lOO; 
y*y0 ( : ) ; 
f*y*zeros (1,6); 
tout»t; 
yout*y . ' ; 

tau*tol*max  (nor(a(y, '  inf' ) ,  1)  ; 

%  Main  loop 

while (t>tf inal) 4 (h>«hmin) 
if  t-h>tfinal  ,h«t-tfinal;  end 
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%  Compute  the  slopes 

temp>feval  (FunFcn,  leal,  tint  1,  alprl,t,  y)  ; 
f  (:,l)>temp(:) ; 
for  j»l:5 

temp»feval (FunFcn, kal, zintl, alprl, t-alpha ( j) *h, y-h*f *beta  ( 
j)); 

f (:,  j+l)*temp(:) ; 
end 

%  Estimate  the  error  and  the  acceptaple  error 
delta-norm (h*f*gamma(:, 2) , ' inf' ) ; 
taustol*max (norm (y, ' inf' ) , 1) ; 

%  Update  the  solution  only  if  the  error  is  acceptaple 
if  delta<>tau 
t«t-h; 

y«y-h*  f  *gamma ( : , 1 ) ; 
tout»{tout;t J ; 
yout»(yout;y.' J ; 
end 

%  Update  the  step  size 
pow«l/5; 
if  delta-^O.O 

h«min  (hmax,  0.8*h*  (tau/delta)  '^pow) ; 

end 

end 

if  (tfinal<t) 

disp(  'singularity  likely.') 
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%  FUNCTION  "TUBE.M"  GIVES  THE  DERIVATIVE  OF  IMPEDANCE 
function  dzclz=tube  (kal,  zintl,  alprl,  zeta,  z) 
dzdz='(  j.  *kal .  *zintl  .*  (l-<z. /zintl) .  *2)  +2.  *alprl  .*z) ; 

%  FUNCTION  "TUBEl.M"  GIVES  THE  DERIVATIVE  OF  PRESSURE 
function  dpdz«tubel (kal, zintl, zdumy, zeta,pl) 
dpdz» j . *pl . *kal . *zintl . /zdumy; 

%  FUNCTION  "ERRO.M"  FOR  "FMINS".  LE/wST  SQUARS  TECHNIQUE. 

function  rniserr«erro(xO,amplit,  freq) 
mainp*xO  (1,1); 
f0«x0(2,l); 
z«x0(3,l); 

foci* (mamp*f 0) . / (z*f req) ; 

err* ( (foci,  *2) ,/  ( (fO./freq-freq/fO)  .''2+l/z'‘2) ) -aroplit . '‘2 
rmserr*suin(err.‘'2) ; 

%  FUNCTION  -ERR03.M"  FOR  "FMIN".  LEAST  SQUARS  TECHNIQUE 
function  fplus«erro3 (wl, itrgh) 

( z, dens, sspeed, vise] *argo3f (wl, itrgh) ; 

gamma*5/3*  (l-»^eps)  ; 

npr*2/3*  (l-»^eps)  ; 

sqri*  ( (1  + j)  /sqrt  (2) )  *  (1-i-eps) ; 


fofw«sqri*sqrt  (dens'‘3*sspeed''4*npr . . . 

. / (wl . *visc) ) . / (gamma-1) ; 
sl*real (fofw) ; 
s2»real  (z) ; 
s3=imag(fofw) ; 
s4*imag(z) ; 

fplus=  (sl*2-s2''2)  *2+  (s3"2-s4"2)  "2; 


%  FUNCTION  "ARG03F.M".  "USED  BY  ERR03.M" 

function 

(z, dens, sspeed, viscI»argo3  <wl, itrgh) 
w*wl; 

%  establish  some  often  used  constants 

s«=5;  %  number  of  sections 

j*sqrt  (-1) ; 
npr=2/3* (1+eps) ; 
gainma»5/3*  (1+eps) ; 
cp»(2.5*8.3143/ (4.e-3) ) * (1+eps) ; 
ycsolid*.16"  (1+eps)  ; 
sqri« ( (1+ j) /sqrt (2) ) * (1+eps) ; 
teren»zeros (s, 1) ; 
numblay^zeros ( 1 , s ) ; 
lengsec*zeros (1,  s) ; 
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tempR»zeros (1,  s) ; 
tempL*zeros (1,  s) ; 
poros*zeros (1,  s) ; 
ratio*zeros (1,  s)  ; 

%  SET  VALUES 

termin^sum (abs (' free' ) > ; 
ainpres«  (1 .  Ole+5)  *  (l-i-eps)  ; 
dramp« (1 .e-8) ; 


%  section  "1" 
terenl=' opentu' ; 
teren (1,1) «sum (abs (terenl ) ) ; 
numblay (1, 1)»1; 

lengsec  (1,1)*  (65.34e-2)  *  (1-feps) ; 
tempRd,  1)  *293*  (l-t-eps) ; 
tempLd,  1)*293*  (1-^eps) ; 
ratio  (1, 1)  *  (1 .917e-2)  *  (l-t-eps) ; 
poros (1, 1) »1* (1+eps) ; 
ares (1, : )*pi*ratio (1, ; ) * (1+eps) ; 
%  section  "2" 
teren2*' hexch' ; 
teren(l,2)*suni(abs(teren2) ) ; 
numblay (1, 2) *1; 

lengsec (1, 2) « ( . 82e-2) * (l+eps) ; 
tempRd,  2)  *293*  (1+eps)  ; 
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tempRd,  5)  =itrgh*  (1+eps)  ; 
tempL (1, 5) *itrgh* (1+eps) ; 
ratio  (1, 5) » (1 . 917e-2) * (1+eps) ; 
poros (1, 5) *1* (1+eps) ; 
if  s>=7 
%  section  "6" 

teren6='hexch' ; 

terend,  6)*sum(abs(teren6) ) ; 

nutnblay  (1,  6)  =1; 

lengsec  (1, 6) *  d .02e“2) * (l+eps) ; 
tempRd,  6)  *293*  (l+eps)  ; 
tempL (1,6) =293*  (1+eps) ; 
ratio  (1, 6) *  (1 .02e~3) * (1+eps) ; 
poros (1, 6) *. 667* (1+eps) ; 

%  section  ••7" 

teren7*' opentu' ; 

teren (1, 7) *sum(abs (teren7) ) ; 

numblay (1, 7) *1; 

lengsec (1,7)* (5.483e-2) * (1+eps) ; 
tempR(l,7)=293* (1+eps); 
tempL (1, 7) =293* (l+eps) ; 
ratio (1, 7) *1 . 917e-2* (l+eps) ; 
poros (1, 7) *1* (l+eps) ; 
end 

%  START  AT  THE  RIGHT  AND  MOVE  AT  THE  LEFT. 
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nuintot»suin(numblay)  -^1; 

dens(l,nuintot)=ampres*4.0e-3  /  (tempRd,  s)  *8 . 3143) ;  %  density 
vise  (1, numtot )  =1 . 887e-5*  (tempRd,  s)  /273. 15)  ^  ( .  6567)  ; 

%  viscosity 

)cgas  (1/ numtot)  =visc  (1#  numtot)  *cp/npr;  %  termal  conductivity 
sspeedd,  numtot)  =972 . 8*sqrt  (tempRd,  s) /273 . 15) ;  %  speed 

%  isothermal 

typel-' free' ; 
type2=' rigid' ; 

if  termin*»sum(abs (type2) ) 

%  IMPEDANCE  OF  RIGID  TERMINATION. 
fac3 <1,  numtot) =sqrt ( (dens (1, numtot) . . . 

*sspeedd, numtot)  .'‘2)  ./  (w.^visc  (1,  numtot) ) ) ; 
z  (1, numtot)  *d+j)  .*  (dens  d, numtot)  .  *sspeedd, numtot)  . . . 

. *fac3 (1, numtot) *sqrt (npr) ) / (sqrt (2) * (gamma-1) ) ; 

%  IMPEDANCE  OF  FREE  TERMINATION 
else 

z (1, numtot) =free (dens (1, numtot) , vise  d, numtot) , w, ratio  (1, s) . 
..  , s speed (1, numtot) , gamma, npr) ; 

end 

%  IMPEDANCE  TRANSLATE  FOR  THE  OPEN  TUBE  OR  HEAT  EXCHANGER, 
fault*'  Sts'"’  ' ; 
truel*' opentu' ; 
true2*'hexch' ; 
for  )c=s:-l:l 
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jup»numt ot -sum { jup+numblay ( 1 , k : s ) ) ; 
if  teren (1, k) -“sumiabs (fault) ) 
densd,  jup)»ampres*4.0e-3  /  (tempR(l,k)  .*8.3143) ;  %  density 
vise  (1,  jup)  »1 . 887e-5.  *  (tempRd,  k)/273.15)  .''(.6567); 

%  viscosity 

kgasd,  jup) -vised,  jup)  .*cp/npr;  %  termal  conductivity 
sspeedd,  jup) -972. 8.*sqrt  (tempRd,  k) /273. 15) ;  %  speed 

%  isothermal 


%  GET  LAMBDA  AND  LAMBDA (T) 
lambda (1,  jup) -ratio (l,k) . *sqrt  (dens (1,  jup)  .*w. /vised,  jup) ) 
lambdt (1, jup)*sqrt (npr) . *lambda (1, jup) ; 

%  GET  F(L),  F(L(T>>  AND  WAVENUMBERS  FOR  OPEN  TUBE 
PARTS 

if  teren(l,k)«sum(abs(truel) ) 
flam (1, jup) -l~(l+j) ,*sqrt (2) ./lambda (1, jup);  %  f(l) 

flamtd,  jup) -l-d+j)  .*sqrt  (2)  ./lambdt  (1,  jup);  %  f(l(t)) 
fac31-(l+ (gamma-1) /sqrt (npr) ) /sqrt (2) ; 
ka  (1,  jup) » (w.  /sspeedd,  jup) ) . . . 

.*  (1+ (1+j) . *fac31. /lambda (1, jup) ) ; 

else 

%  GET  F(L),  F(L(T))  AND  WAVENUMBERS  FOR  HEAT 
EXCHANGERS  TYPE  "SLIT*. 

sqrmi- (1-j) /sqrt (2) ; 

ar ( 1 , jup) -sqrmi . * lambda ( 1 , jup) ; 
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argumd,  jup)  «exp  ar  (1,  jup) )  ; 
ar  (1,  jup)  *ai  (1,  jup)  /2; 

ctanhd,  jup)  ®  (l-argumd,  jup) )  ./  d+argurad,  jup) )  / 
flam  (1,  jup)  »l“Ctanh  d,  jup) . /ar  (1,  jup) ;  %  fd) 

ar (1,  jup) »sqrmi*lambdt (1, jup) ; 
argumd,  jup)*exp(-ar  (1,  jup) )  ; 
ard,  jup)«ard,  jup)/2; 

ct anh ( 1 , jup) * ( 1-argum ( 1 , jup) ) . / ( 1 +argum ( 1 , jup) ) ; 
flamt (1, jup) =l-ctanh d, jup) . /ar (1, jup) ;  %  f(l(t)) 
low»w/sspeed ( 1 ,  jup) ; 

kail,  jup) »low . *sqrt ( (gamma- (gamma-1 ) . . . 

.  *flamt  (1,  jup) )  .  /flamd,  jup) )  ; 

end 

%  TRANSLATION  THEOREM 
comden  (1,  jup) «dens  (1,  jup) . /  flamd,  jup) ; 
zint  (1,  jup)  -cornden  (1,  jup) .  *w./  (Jca(l,  jup) .  *poros  (l,)c) ) ; 
dsub (1, k) "lengsec (1, k) . /numblay (1,  k) ; 
ar (1, jup) >ka (1, jup) . *dsub(l,  k) ; 
snd,  jup)*sin(ar  (1,  jup) )  ; 
csd,  jup)«cos(ard,  jup)); 
ct  (1,  jup)»cs  (1,  jup)  ./snd,  jup)  ; 
fac3  (1,  jup)  «zint  (1,  jup) .  /z  (1, 1  + jup) ; 
z (1, jup) -zint (1, jup) (ct (1, jup) . . . 

- j. *fac3 (1, jup) ) ./ (fac3(l, jup) .*ct (1, jup)-j) ; 

else 
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%  STACK  WITH  LINEAR  DISTRIBUTION  BETWEEN  THE  ENDS  OF  THE 
STACK. 

t02=(tempR(l,k)“tempL(l,k) ) /lengsec(l,lc) ;  %  approx 

dz/dt 

ns=0; 

for  la=k+numblay (1, k) -1 :-l : k 
ns®ns+l; 

tens“tempR(l,k)-(terapR(l,k)-tempL(l,k) )  *ns/nuinblay  (1, k) ; 
tnsml=tempR(l,k)-(tempR(l,k)-tempL(l,k) ) * (ns-1) /numblay (1, k) 

i 

tave (1, la) = (tens+tnsml) /2; 

dsub (1, k) =lengsec (1, k) /numblay k) ; 

dens (1, la) =ampres*4 . Oe-3  /  (tave (1, la) *8. 3143) ;  %  density 

vise (1, la) =1 . 887e-5*  (tave  (1, la) /273. 15) " ( . 6567) ;  % 

viscosity 

kgas (1, la) «visc (1, jup) *cp/npr;  %  termal  conductivity 
sspeed(l,la)»972.8*sqrt  (taved,  la) /273.15) ;  %  speed 

%  isothermal 

%  GET  LAMBDA  AND  LAMBDA (T) 
lambda (1, la) >ratio (l,k) . *sqrt (dens (1, la) . *w. /vise (1, la) ) ; 
lambdt (1, la) =sqrt (npr) ^lambda (1,  la) ; 
ar (1, la) ^sqrmi*lambda (1, la) ; 
argum (1, la) =exp (-ar (1, la) ) ; 
ar (1, la) ^ar (1, la) /2; 

ctanh  (1,  la) » (1-argum  (1,  la) ) .  /  (l-»-argum  (1,  la) )  ; 
flamd,  la)  =l-ctanh  (1,  la) . /ar  (1,  la) ;  %  fd) 
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ar  (1,  la)  «sqrmi*lambcit  (1,  la) ; 
argum (1, la) =exp (-ar  <1, la) ) ; 
ar  (1, la) *ar (1, la) /2; 

ctanh  (1,  la)  =  (l-argumd,  la) )  .  /  (l+argum  (1,  la) )  ; 
flamt(l,la)»l-ctanh(l,la)  ./ar(l,la);  %  f(l(t)) 

1 ow>w/s speed (1, la) ; 

)ca  (1,  la)  *low.  *sqrt  ( (ganuna-  (gaimna-1)  . . . 

.  *flamt  (1,  la) ) .  /flamd,  la) ) ; 

zint (1, la) =dens (1, la) . *w _ 

.  /  (poros  (1,  )c)  .  '►flarod,  la)  .  *Jca  (1,  la) )  ; 

alprimd,  la)  «toz*  (flamt  (1,  la) _ 

.  /flam  (1,  la) -1)  .  /  (2*tave  (1,  la)  *  d-npr) )  ; 
)cal-)ca  (1,  la)  .  * ; 
zintl*zint (1, la) 
alprl=Balprimd,  la) .  * ; 

zetaO=  dengsec  d,k)-  (ns-1)  *lengsec  d,  Jc)  /numblay  (1, )c) ) ; 
zetafin=  (lengsec  d,  )c)-  (ns)  *lengsec  (1,  k)  /numblay  d,  k) )  ; 
zO=z (1, 1+la) . ' ; 
tol»l .e-4; 

(zeta, zout ]*trode45 ('tube' , zetaO, zetafin,  zO,  tol,  kal, zintl, al 
prl)  ; 

n=length (zout) ; 
z (1, la) «zout (n, 1) ; 
end 
end 
end 
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[•] 


z=z (1,1); 
dens=ciens  (1,1); 
sspeed*sspeecl  (1,1); 
visc=visc (1, 1) ; 
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